/*
 * @file bench_eigen_mtimes.cpp
 *
 * Benchmarking of Eigen matrix multiplication
 *
 * @author Dahua Lin
 */


#include <gmlpp/base/eigen_base.h>
#include "bench_base.h"

using namespace gml;

class MatXdMtimesBench : public BenchTask<MatXdMtimesBench>
{
public:
	MatXdMtimesBench(int n)
	: m_dim(n)
	, A(Eigen::MatrixXd::Random(n, n))
	, B(Eigen::MatrixXd::Random(n, n))
	, C(Eigen::MatrixXd::Zero(n, n))
	{ }

	GML_ENSURE_INLINE void run()
	{
		C.noalias() = A * B;
	}

private:
	int m_dim;
	Eigen::MatrixXd A;
	Eigen::MatrixXd B;
	Eigen::MatrixXd C;
};



template<int Dim>
class MatFdMtimesBench : public BenchTask<MatFdMtimesBench<Dim> >
{
public:
	typedef Eigen::Matrix<double, Dim, Dim> MatT;

	MatFdMtimesBench()
	: A(MatT::Random())
	, B(MatT::Random())
	, C(MatT::Zero())
	{ }

	GML_ENSURE_INLINE void run()
	{
		C.noalias() = A * B;
	}

private:
	MatT A;
	MatT B;
	MatT C;
};




static double GLOBAL_MKL_BASE_GFLOPS = 1;

inline void assess_gflops()
{
	MatXdMtimesBench tsk(1000);
	BenchStats ss = bench_task(tsk, 1, 1);
	GLOBAL_MKL_BASE_GFLOPS = 2 * ss.rate();
}


template<class T>
void do_bench(BenchTask<T>& task, size_t n, size_t ntimes)
{
	BenchStats ss = bench_task(task, ntimes, 1);

	double nd = double(n);
	double gflops = 2 * (nd * nd * nd) * 1.0e-9 * ss.rate();

	std::printf("    n = %4lu:   %8.4f GFLOPS  | rate = %7.2e/s  elapsed = %5.2f s\n",
			n, gflops, ss.rate(), ss.elapsed_seconds);
}


size_t calc_ntimes_mkl(int n, double gbase)
{
	double base_gflops = gbase;
	if (n <= 64) base_gflops *= 0.75;
	if (n <= 32) base_gflops *= 0.75;
	if (n <= 16) base_gflops *= 0.75;
	if (n <= 8) base_gflops *= 0.5;
	if (n <= 6) base_gflops *= 0.5;
	if (n <= 4) base_gflops *= 0.5;
	if (n <= 3) base_gflops *= 0.3;
	if (n <= 2) base_gflops *= 0.3;

	double nd = double(n);
	double gpt = (2 * nd * nd * nd) * 1.0e-9;

	return size_t(base_gflops / gpt * 0.5) + 1;
}


size_t calc_ntimes_sta(int n, double gbase)
{
	double base_gflops = gbase * 0.5;

	double nd = double(n);
	double gpt = (2 * nd * nd * nd) * 1.0e-9;

	return size_t(base_gflops / gpt * 0.5) + 1;
}




void do_bench_dyn(int n)
{
	size_t ntimes = calc_ntimes_mkl(n, GLOBAL_MKL_BASE_GFLOPS);
	MatXdMtimesBench task(n);
	do_bench(task, n, ntimes);
}

template<int Dim>
void do_bench_sta()
{
	size_t ntimes = calc_ntimes_sta(Dim, GLOBAL_MKL_BASE_GFLOPS);

	MatFdMtimesBench<Dim> task;
	do_bench(task, Dim, ntimes);
}


int main(int argc, char *argv[])
{
	std::printf("Matrix Multiplication Benchmark\n");
	std::printf("*************************************\n");

	std::printf("Performance calibration ...\n");
	assess_gflops();
	std::printf("\n");

	std::printf("Static size:\n");
	std::printf("---------------------------\n");

	do_bench_sta<4>();
	do_bench_sta<6>();
	do_bench_sta<8>();
	do_bench_sta<12>();
	do_bench_sta<16>();

	std::printf("\n");

	std::printf("Dynamic size:\n");
	std::printf("---------------------------\n");

	int mkl_test_dims[] = {
			2, 3, 4, 5, 6, 8, 10,
			12, 16, 20, 25, 30, 40, 60, 80,
			100, 150, 200, 300, 400, 600, 800,
			1000, 1200, 1600 };

	size_t num_mkl_test_dims = sizeof(mkl_test_dims) / sizeof(int);

	for (size_t k = 0; k < num_mkl_test_dims; ++k)
	{
		int dim = mkl_test_dims[k];
		do_bench_dyn(dim);
	}

	std::printf("\n");

	return 0;
}




